{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "google",
   "metadata": {},
   "source": [
    "##### Copyright 2023 Google LLC."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "apache",
   "metadata": {},
   "source": [
    "Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "you may not use this file except in compliance with the License.\n",
    "You may obtain a copy of the License at\n",
    "\n",
    "    http://www.apache.org/licenses/LICENSE-2.0\n",
    "\n",
    "Unless required by applicable law or agreed to in writing, software\n",
    "distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "See the License for the specific language governing permissions and\n",
    "limitations under the License.\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "basename",
   "metadata": {},
   "source": [
    "# rcpsp_sat"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "link",
   "metadata": {},
   "source": [
    "<table align=\"left\">\n",
    "<td>\n",
    "<a href=\"https://colab.research.google.com/github/google/or-tools/blob/main/examples/notebook/examples/rcpsp_sat.ipynb\"><img src=\"https://raw.githubusercontent.com/google/or-tools/main/tools/colab_32px.png\"/>Run in Google Colab</a>\n",
    "</td>\n",
    "<td>\n",
    "<a href=\"https://github.com/google/or-tools/blob/main/examples/python/rcpsp_sat.py\"><img src=\"https://raw.githubusercontent.com/google/or-tools/main/tools/github_32px.png\"/>View source on GitHub</a>\n",
    "</td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "doc",
   "metadata": {},
   "source": [
    "First, you must install [ortools](https://pypi.org/project/ortools/) package in this colab."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "install",
   "metadata": {},
   "outputs": [],
   "source": [
    "%pip install ortools"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "description",
   "metadata": {},
   "source": [
    "\n",
    "Sat based solver for the RCPSP problems (see rcpsp.proto).\n",
    "\n",
    "Introduction to the problem:\n",
    "   https://www.projectmanagement.ugent.be/research/project_scheduling/rcpsp\n",
    "\n",
    "Data use in flags:\n",
    "  http://www.om-db.wi.tum.de/psplib/data.html\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "code",
   "metadata": {},
   "outputs": [],
   "source": [
    "import collections\n",
    "\n",
    "from ortools.sat.colab import flags\n",
    "from google.protobuf import text_format\n",
    "from ortools.sat.python import cp_model\n",
    "from ortools.scheduling import rcpsp_pb2\n",
    "from ortools.scheduling.python import rcpsp\n",
    "\n",
    "_INPUT = flags.define_string(\"input\", \"\", \"Input file to parse and solve.\")\n",
    "_OUTPUT_PROTO = flags.define_string(\n",
    "    \"output_proto\", \"\", \"Output file to write the cp_model proto to.\"\n",
    ")\n",
    "_PARAMS = flags.define_string(\"params\", \"\", \"Sat solver parameters.\")\n",
    "_USE_INTERVAL_MAKESPAN = flags.define_bool(\n",
    "    \"use_interval_makespan\",\n",
    "    True,\n",
    "    \"Whether we encode the makespan using an interval or not.\",\n",
    ")\n",
    "_HORIZON = flags.define_integer(\"horizon\", -1, \"Force horizon.\")\n",
    "_ADD_REDUNDANT_ENERGETIC_CONSTRAINTS = flags.define_bool(\n",
    "    \"add_redundant_energetic_constraints\",\n",
    "    False,\n",
    "    \"Add redundant energetic constraints on the pairs of tasks extracted from\"\n",
    "    + \" precedence graph.\",\n",
    ")\n",
    "_DELAY_TIME_LIMIT = flags.define_float(\n",
    "    \"delay_time_limit\",\n",
    "    20.0,\n",
    "    \"Time limit when computing min delay between tasks.\"\n",
    "    + \" A non-positive time limit disable min delays computation.\",\n",
    ")\n",
    "_PREEMPTIVE_LB_TIME_LIMIT = flags.define_float(\n",
    "    \"preemptive_lb_time_limit\",\n",
    "    0.0,\n",
    "    \"Time limit when computing a preemptive schedule lower bound.\"\n",
    "    + \" A non-positive time limit disable this computation.\",\n",
    ")\n",
    "\n",
    "\n",
    "def PrintProblemStatistics(problem):\n",
    "    \"\"\"Display various statistics on the problem.\"\"\"\n",
    "\n",
    "    # Determine problem type.\n",
    "    problem_type = (\n",
    "        \"Resource Investment Problem\" if problem.is_resource_investment else \"RCPSP\"\n",
    "    )\n",
    "\n",
    "    num_resources = len(problem.resources)\n",
    "    num_tasks = len(problem.tasks) - 2  # 2 sentinels.\n",
    "    tasks_with_alternatives = 0\n",
    "    variable_duration_tasks = 0\n",
    "    tasks_with_delay = 0\n",
    "\n",
    "    for task in problem.tasks:\n",
    "        if len(task.recipes) > 1:\n",
    "            tasks_with_alternatives += 1\n",
    "            duration_0 = task.recipes[0].duration\n",
    "            for recipe in task.recipes:\n",
    "                if recipe.duration != duration_0:\n",
    "                    variable_duration_tasks += 1\n",
    "                    break\n",
    "        if task.successor_delays:\n",
    "            tasks_with_delay += 1\n",
    "\n",
    "    if problem.is_rcpsp_max:\n",
    "        problem_type += \"/Max delay\"\n",
    "    # We print 2 less tasks as these are sentinel tasks that are not counted in\n",
    "    # the description of the rcpsp models.\n",
    "    if problem.is_consumer_producer:\n",
    "        print(f\"Solving {problem_type} with:\")\n",
    "        print(f\"  - {num_resources} reservoir resources\")\n",
    "        print(f\"  - {num_tasks} tasks\")\n",
    "    else:\n",
    "        print(f\"Solving {problem_type} with:\")\n",
    "        print(f\"  - {num_resources} renewable resources\")\n",
    "        print(f\"  - {num_tasks} tasks\")\n",
    "        if tasks_with_alternatives:\n",
    "            print(f\"    - {tasks_with_alternatives} tasks with alternative resources\")\n",
    "        if variable_duration_tasks:\n",
    "            print(f\"    - {variable_duration_tasks} tasks with variable durations\")\n",
    "        if tasks_with_delay:\n",
    "            print(f\"    - {tasks_with_delay} tasks with successor delays\")\n",
    "\n",
    "\n",
    "def AnalyseDependencyGraph(problem):\n",
    "    \"\"\"Analyses the dependency graph to improve the model.\n",
    "\n",
    "    Args:\n",
    "      problem: the protobuf of the problem to solve.\n",
    "\n",
    "    Returns:\n",
    "      a list of (task1, task2, in_between_tasks) with task2 and indirect successor\n",
    "      of task1, and in_between_tasks being the list of all tasks after task1 and\n",
    "      before task2.\n",
    "    \"\"\"\n",
    "\n",
    "    num_nodes = len(problem.tasks)\n",
    "    print(f\"Analysing the dependency graph over {num_nodes} nodes\")\n",
    "\n",
    "    ins = collections.defaultdict(list)\n",
    "    outs = collections.defaultdict(list)\n",
    "    after = collections.defaultdict(set)\n",
    "    before = collections.defaultdict(set)\n",
    "\n",
    "    # Build the transitive closure of the precedences.\n",
    "    # This algorithm has the wrong complexity (n^4), but is OK for the psplib\n",
    "    # as the biggest example has 120 nodes.\n",
    "    for n in range(num_nodes):\n",
    "        for s in problem.tasks[n].successors:\n",
    "            ins[s].append(n)\n",
    "            outs[n].append(s)\n",
    "\n",
    "            for a in list(after[s]) + [s]:\n",
    "                for b in list(before[n]) + [n]:\n",
    "                    after[b].add(a)\n",
    "                    before[a].add(b)\n",
    "\n",
    "    # Search for pair of tasks, containing at least two parallel branch between\n",
    "    # them in the precedence graph.\n",
    "    num_candidates = 0\n",
    "    result = []\n",
    "    for source, start_outs in outs.items():\n",
    "        if len(start_outs) <= 1:\n",
    "            # Starting with the unique successor of source will be as good.\n",
    "            continue\n",
    "        for sink, end_ins in ins.items():\n",
    "            if len(end_ins) <= 1:\n",
    "                # Ending with the unique predecessor of sink will be as good.\n",
    "                continue\n",
    "            if sink == source:\n",
    "                continue\n",
    "            if sink not in after[source]:\n",
    "                continue\n",
    "\n",
    "            num_active_outgoing_branches = 0\n",
    "            num_active_incoming_branches = 0\n",
    "            for succ in outs[source]:\n",
    "                if sink in after[succ]:\n",
    "                    num_active_outgoing_branches += 1\n",
    "            for pred in ins[sink]:\n",
    "                if source in before[pred]:\n",
    "                    num_active_incoming_branches += 1\n",
    "\n",
    "            if num_active_outgoing_branches <= 1 or num_active_incoming_branches <= 1:\n",
    "                continue\n",
    "\n",
    "            common = after[source].intersection(before[sink])\n",
    "            if len(common) <= 1:\n",
    "                continue\n",
    "            num_candidates += 1\n",
    "            result.append((source, sink, common))\n",
    "\n",
    "    # Sort entries lexicographically by (len(common), source, sink)\n",
    "    def Price(entry):\n",
    "        return num_nodes * num_nodes * len(entry[2]) + num_nodes * entry[0] + entry[1]\n",
    "\n",
    "    result.sort(key=Price)\n",
    "    print(f\"  - created {len(result)} pairs of nodes to examine\", flush=True)\n",
    "    return result, after\n",
    "\n",
    "\n",
    "def SolveRcpsp(\n",
    "    problem,\n",
    "    proto_file,\n",
    "    params,\n",
    "    active_tasks,\n",
    "    source,\n",
    "    sink,\n",
    "    intervals_of_tasks,\n",
    "    delays,\n",
    "    in_main_solve=False,\n",
    "    initial_solution=None,\n",
    "    lower_bound=0,\n",
    "):\n",
    "    \"\"\"Parse and solve a given RCPSP problem in proto format.\n",
    "\n",
    "    The model will only look at the tasks {source} + {sink} + active_tasks, and\n",
    "    ignore all others.\n",
    "\n",
    "    Args:\n",
    "      problem: the description of the model to solve in protobuf format\n",
    "      proto_file: the name of the file to export the CpModel proto to.\n",
    "      params: the string representation of the parameters to pass to the sat\n",
    "        solver.\n",
    "      active_tasks: the set of active tasks to consider.\n",
    "      source: the source task in the graph. Its end will be forced to 0.\n",
    "      sink: the sink task of the graph. Its start is the makespan of the problem.\n",
    "      intervals_of_tasks: a heuristic lists of (task1, task2, tasks) used to add\n",
    "        redundant energetic equations to the model.\n",
    "      delays: a list of (task1, task2, min_delays) used to add extended precedence\n",
    "        constraints (start(task2) >= end(task1) + min_delay).\n",
    "      in_main_solve: indicates if this is the main solve procedure.\n",
    "      initial_solution: A valid assignment used to hint the search.\n",
    "      lower_bound: A valid lower bound of the makespan objective.\n",
    "\n",
    "    Returns:\n",
    "      (lower_bound of the objective, best solution found, assignment)\n",
    "    \"\"\"\n",
    "    # Create the model.\n",
    "    model = cp_model.CpModel()\n",
    "    model.SetName(problem.name)\n",
    "\n",
    "    num_resources = len(problem.resources)\n",
    "\n",
    "    all_active_tasks = list(active_tasks)\n",
    "    all_active_tasks.sort()\n",
    "    all_resources = range(num_resources)\n",
    "\n",
    "    horizon = problem.deadline if problem.deadline != -1 else problem.horizon\n",
    "    if _HORIZON.value > 0:\n",
    "        horizon = _HORIZON.value\n",
    "    elif delays and in_main_solve and (source, sink) in delays:\n",
    "        horizon = delays[(source, sink)][1]\n",
    "    elif horizon == -1:  # Naive computation.\n",
    "        horizon = sum(max(r.duration for r in t.recipes) for t in problem.tasks)\n",
    "        if problem.is_rcpsp_max:\n",
    "            for t in problem.tasks:\n",
    "                for sd in t.successor_delays:\n",
    "                    for rd in sd.recipe_delays:\n",
    "                        for d in rd.min_delays:\n",
    "                            horizon += abs(d)\n",
    "    if in_main_solve:\n",
    "        print(f\"Horizon = {horizon}\", flush=True)\n",
    "\n",
    "    # Containers.\n",
    "    task_starts = {}\n",
    "    task_ends = {}\n",
    "    task_durations = {}\n",
    "    task_intervals = {}\n",
    "    task_resource_to_energy = {}\n",
    "    task_to_resource_demands = collections.defaultdict(list)\n",
    "\n",
    "    task_to_presence_literals = collections.defaultdict(list)\n",
    "    task_to_recipe_durations = collections.defaultdict(list)\n",
    "    task_resource_to_fixed_demands = collections.defaultdict(dict)\n",
    "    task_resource_to_max_energy = collections.defaultdict(int)\n",
    "\n",
    "    resource_to_sum_of_demand_max = collections.defaultdict(int)\n",
    "\n",
    "    # Create task variables.\n",
    "    for t in all_active_tasks:\n",
    "        task = problem.tasks[t]\n",
    "        num_recipes = len(task.recipes)\n",
    "        all_recipes = range(num_recipes)\n",
    "\n",
    "        start_var = model.NewIntVar(0, horizon, f\"start_of_task_{t}\")\n",
    "        end_var = model.NewIntVar(0, horizon, f\"end_of_task_{t}\")\n",
    "\n",
    "        literals = []\n",
    "        if num_recipes > 1:\n",
    "            # Create one literal per recipe.\n",
    "            literals = [model.NewBoolVar(f\"is_present_{t}_{r}\") for r in all_recipes]\n",
    "\n",
    "            # Exactly one recipe must be performed.\n",
    "            model.AddExactlyOne(literals)\n",
    "\n",
    "        else:\n",
    "            literals = [1]\n",
    "\n",
    "        # Temporary data structure to fill in 0 demands.\n",
    "        demand_matrix = collections.defaultdict(int)\n",
    "\n",
    "        # Scan recipes and build the demand matrix and the vector of durations.\n",
    "        for recipe_index, recipe in enumerate(task.recipes):\n",
    "            task_to_recipe_durations[t].append(recipe.duration)\n",
    "            for demand, resource in zip(recipe.demands, recipe.resources):\n",
    "                demand_matrix[(resource, recipe_index)] = demand\n",
    "\n",
    "        # Create the duration variable from the accumulated durations.\n",
    "        duration_var = model.NewIntVarFromDomain(\n",
    "            cp_model.Domain.FromValues(task_to_recipe_durations[t]),\n",
    "            f\"duration_of_task_{t}\",\n",
    "        )\n",
    "\n",
    "        # Link the recipe literals and the duration_var.\n",
    "        for r in range(num_recipes):\n",
    "            model.Add(duration_var == task_to_recipe_durations[t][r]).OnlyEnforceIf(\n",
    "                literals[r]\n",
    "            )\n",
    "\n",
    "        # Create the interval of the task.\n",
    "        task_interval = model.NewIntervalVar(\n",
    "            start_var, duration_var, end_var, f\"task_interval_{t}\"\n",
    "        )\n",
    "\n",
    "        # Store task variables.\n",
    "        task_starts[t] = start_var\n",
    "        task_ends[t] = end_var\n",
    "        task_durations[t] = duration_var\n",
    "        task_intervals[t] = task_interval\n",
    "        task_to_presence_literals[t] = literals\n",
    "\n",
    "        # Create the demand variable of the task for each resource.\n",
    "        for res in all_resources:\n",
    "            demands = [demand_matrix[(res, recipe)] for recipe in all_recipes]\n",
    "            task_resource_to_fixed_demands[(t, res)] = demands\n",
    "            demand_var = model.NewIntVarFromDomain(\n",
    "                cp_model.Domain.FromValues(demands), f\"demand_{t}_{res}\"\n",
    "            )\n",
    "            task_to_resource_demands[t].append(demand_var)\n",
    "\n",
    "            # Link the recipe literals and the demand_var.\n",
    "            for r in all_recipes:\n",
    "                model.Add(demand_var == demand_matrix[(res, r)]).OnlyEnforceIf(\n",
    "                    literals[r]\n",
    "                )\n",
    "\n",
    "            resource_to_sum_of_demand_max[res] += max(demands)\n",
    "\n",
    "        # Create the energy expression for (task, resource):\n",
    "        for res in all_resources:\n",
    "            task_resource_to_energy[(t, res)] = sum(\n",
    "                literals[r]\n",
    "                * task_to_recipe_durations[t][r]\n",
    "                * task_resource_to_fixed_demands[(t, res)][r]\n",
    "                for r in all_recipes\n",
    "            )\n",
    "            task_resource_to_max_energy[(t, res)] = max(\n",
    "                task_to_recipe_durations[t][r]\n",
    "                * task_resource_to_fixed_demands[(t, res)][r]\n",
    "                for r in all_recipes\n",
    "            )\n",
    "\n",
    "    # Create makespan variable\n",
    "    makespan = model.NewIntVar(lower_bound, horizon, \"makespan\")\n",
    "    makespan_size = model.NewIntVar(1, horizon, \"interval_makespan_size\")\n",
    "    interval_makespan = model.NewIntervalVar(\n",
    "        makespan, makespan_size, model.NewConstant(horizon + 1), \"interval_makespan\"\n",
    "    )\n",
    "\n",
    "    # Add precedences.\n",
    "    if problem.is_rcpsp_max:\n",
    "        # In RCPSP/Max problem, precedences are given and max delay (possible\n",
    "        # negative) between the starts of two tasks.\n",
    "        for task_id in all_active_tasks:\n",
    "            task = problem.tasks[task_id]\n",
    "            num_modes = len(task.recipes)\n",
    "\n",
    "            for successor_index, next_id in enumerate(task.successors):\n",
    "                delay_matrix = task.successor_delays[successor_index]\n",
    "                num_next_modes = len(problem.tasks[next_id].recipes)\n",
    "                for m1 in range(num_modes):\n",
    "                    s1 = task_starts[task_id]\n",
    "                    p1 = task_to_presence_literals[task_id][m1]\n",
    "                    if next_id == sink:\n",
    "                        delay = delay_matrix.recipe_delays[m1].min_delays[0]\n",
    "                        model.Add(s1 + delay <= makespan).OnlyEnforceIf(p1)\n",
    "                    else:\n",
    "                        for m2 in range(num_next_modes):\n",
    "                            delay = delay_matrix.recipe_delays[m1].min_delays[m2]\n",
    "                            s2 = task_starts[next_id]\n",
    "                            p2 = task_to_presence_literals[next_id][m2]\n",
    "                            model.Add(s1 + delay <= s2).OnlyEnforceIf([p1, p2])\n",
    "    else:\n",
    "        # Normal dependencies (task ends before the start of successors).\n",
    "        for t in all_active_tasks:\n",
    "            for n in problem.tasks[t].successors:\n",
    "                if n == sink:\n",
    "                    model.Add(task_ends[t] <= makespan)\n",
    "                elif n in active_tasks:\n",
    "                    model.Add(task_ends[t] <= task_starts[n])\n",
    "\n",
    "    # Containers for resource investment problems.\n",
    "    capacities = []  # Capacity variables for all resources.\n",
    "    max_cost = 0  # Upper bound on the investment cost.\n",
    "\n",
    "    # Create resources.\n",
    "    for res in all_resources:\n",
    "        resource = problem.resources[res]\n",
    "        c = resource.max_capacity\n",
    "        if c == -1:\n",
    "            print(f\"No capacity: {resource}\")\n",
    "            c = resource_to_sum_of_demand_max[res]\n",
    "\n",
    "        # RIP problems have only renewable resources, and no makespan.\n",
    "        if problem.is_resource_investment or resource.renewable:\n",
    "            intervals = [task_intervals[t] for t in all_active_tasks]\n",
    "            demands = [task_to_resource_demands[t][res] for t in all_active_tasks]\n",
    "\n",
    "            if problem.is_resource_investment:\n",
    "                capacity = model.NewIntVar(0, c, f\"capacity_of_{res}\")\n",
    "                model.AddCumulative(intervals, demands, capacity)\n",
    "                capacities.append(capacity)\n",
    "                max_cost += c * resource.unit_cost\n",
    "            else:  # Standard renewable resource.\n",
    "                if _USE_INTERVAL_MAKESPAN.value:\n",
    "                    intervals.append(interval_makespan)\n",
    "                    demands.append(c)\n",
    "\n",
    "                model.AddCumulative(intervals, demands, c)\n",
    "        else:  # Non empty non renewable resource. (single mode only)\n",
    "            if problem.is_consumer_producer:\n",
    "                reservoir_starts = []\n",
    "                reservoir_demands = []\n",
    "                for t in all_active_tasks:\n",
    "                    if task_resource_to_fixed_demands[(t, res)][0]:\n",
    "                        reservoir_starts.append(task_starts[t])\n",
    "                        reservoir_demands.append(\n",
    "                            task_resource_to_fixed_demands[(t, res)][0]\n",
    "                        )\n",
    "                model.AddReservoirConstraint(\n",
    "                    reservoir_starts,\n",
    "                    reservoir_demands,\n",
    "                    resource.min_capacity,\n",
    "                    resource.max_capacity,\n",
    "                )\n",
    "            else:  # No producer-consumer. We just sum the demands.\n",
    "                model.Add(\n",
    "                    cp_model.LinearExpr.Sum(\n",
    "                        [task_to_resource_demands[t][res] for t in all_active_tasks]\n",
    "                    )\n",
    "                    <= c\n",
    "                )\n",
    "\n",
    "    # Objective.\n",
    "    if problem.is_resource_investment:\n",
    "        objective = model.NewIntVar(0, max_cost, \"capacity_costs\")\n",
    "        model.Add(\n",
    "            objective\n",
    "            == sum(\n",
    "                problem.resources[i].unit_cost * capacities[i]\n",
    "                for i in range(len(capacities))\n",
    "            )\n",
    "        )\n",
    "    else:\n",
    "        objective = makespan\n",
    "\n",
    "    model.Minimize(objective)\n",
    "\n",
    "    # Add min delay constraints.\n",
    "    if delays is not None:\n",
    "        for (local_start, local_end), (min_delay, _) in delays.items():\n",
    "            if local_start == source and local_end in active_tasks:\n",
    "                model.Add(task_starts[local_end] >= min_delay)\n",
    "            elif local_start in active_tasks and local_end == sink:\n",
    "                model.Add(makespan >= task_ends[local_start] + min_delay)\n",
    "            elif local_start in active_tasks and local_end in active_tasks:\n",
    "                model.Add(task_starts[local_end] >= task_ends[local_start] + min_delay)\n",
    "\n",
    "    problem_is_single_mode = True\n",
    "    for t in all_active_tasks:\n",
    "        if len(task_to_presence_literals[t]) > 1:\n",
    "            problem_is_single_mode = False\n",
    "            break\n",
    "\n",
    "    # Add sentinels.\n",
    "    task_starts[source] = 0\n",
    "    task_ends[source] = 0\n",
    "    task_to_presence_literals[0].append(True)\n",
    "    task_starts[sink] = makespan\n",
    "    task_to_presence_literals[sink].append(True)\n",
    "\n",
    "    # For multi-mode problems, add a redundant energetic constraint:\n",
    "    # for every (start, end, in_between_tasks) extracted from the precedence\n",
    "    # graph, it add the energetic relaxation:\n",
    "    #   (start_var('end') - end_var('start')) * capacity_max >=\n",
    "    #        sum of linearized energies of all tasks from 'in_between_tasks'\n",
    "    if (\n",
    "        not problem.is_resource_investment\n",
    "        and not problem.is_consumer_producer\n",
    "        and _ADD_REDUNDANT_ENERGETIC_CONSTRAINTS.value\n",
    "        and in_main_solve\n",
    "        and not problem_is_single_mode\n",
    "    ):\n",
    "        added_constraints = 0\n",
    "        ignored_constraits = 0\n",
    "        for local_start, local_end, common in intervals_of_tasks:\n",
    "            for res in all_resources:\n",
    "                resource = problem.resources[res]\n",
    "                if not resource.renewable:\n",
    "                    continue\n",
    "                c = resource.max_capacity\n",
    "                if delays and (local_start, local_end) in delays:\n",
    "                    min_delay, _ = delays[local_start, local_end]\n",
    "                    sum_of_max_energies = sum(\n",
    "                        task_resource_to_max_energy[(t, res)] for t in common\n",
    "                    )\n",
    "                    if sum_of_max_energies <= c * min_delay:\n",
    "                        ignored_constraits += 1\n",
    "                        continue\n",
    "                model.Add(\n",
    "                    c * (task_starts[local_end] - task_ends[local_start])\n",
    "                    >= sum(task_resource_to_energy[(t, res)] for t in common)\n",
    "                )\n",
    "                added_constraints += 1\n",
    "        print(\n",
    "            f\"Added {added_constraints} redundant energetic constraints, and \"\n",
    "            + f\"ignored {ignored_constraits} constraints.\",\n",
    "            flush=True,\n",
    "        )\n",
    "\n",
    "    # Add solution hint.\n",
    "    if initial_solution:\n",
    "        for t in all_active_tasks:\n",
    "            model.AddHint(task_starts[t], initial_solution.start_of_task[t])\n",
    "            if len(task_to_presence_literals[t]) > 1:\n",
    "                selected = initial_solution.selected_recipe_of_task[t]\n",
    "                model.AddHint(task_to_presence_literals[t][selected], 1)\n",
    "\n",
    "    # Write model to file.\n",
    "    if proto_file:\n",
    "        print(f\"Writing proto to{proto_file}\")\n",
    "        model.ExportToFile(proto_file)\n",
    "\n",
    "    # Solve model.\n",
    "    solver = cp_model.CpSolver()\n",
    "\n",
    "    # Parse user specified parameters.\n",
    "    if params:\n",
    "        text_format.Parse(params, solver.parameters)\n",
    "\n",
    "    # Favor objective_shaving_search over objective_lb_search.\n",
    "    if solver.parameters.num_workers >= 16 and solver.parameters.num_workers < 24:\n",
    "        solver.parameters.ignore_subsolvers.append(\"objective_lb_search\")\n",
    "        solver.parameters.extra_subsolvers.append(\"objective_shaving_search\")\n",
    "\n",
    "    # Experimental: Specify the fact that the objective is a makespan\n",
    "    solver.parameters.push_all_tasks_toward_start = True\n",
    "\n",
    "    # Enable logging in the main solve.\n",
    "\n",
    "    if in_main_solve:\n",
    "        solver.parameters.log_search_progress = True\n",
    "    #\n",
    "    status = solver.Solve(model)\n",
    "    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:\n",
    "        assignment = rcpsp_pb2.RcpspAssignment()\n",
    "        for t, _ in enumerate(problem.tasks):\n",
    "            if t in task_starts:\n",
    "                assignment.start_of_task.append(solver.Value(task_starts[t]))\n",
    "                for r, recipe_literal in enumerate(task_to_presence_literals[t]):\n",
    "                    if solver.BooleanValue(recipe_literal):\n",
    "                        assignment.selected_recipe_of_task.append(r)\n",
    "                        break\n",
    "            else:  # t is not an active task.\n",
    "                assignment.start_of_task.append(0)\n",
    "                assignment.selected_recipe_of_task.append(0)\n",
    "        return (\n",
    "            int(solver.BestObjectiveBound()),\n",
    "            int(solver.ObjectiveValue()),\n",
    "            assignment,\n",
    "        )\n",
    "    return -1, -1, None\n",
    "\n",
    "\n",
    "def ComputeDelaysBetweenNodes(problem, task_intervals):\n",
    "    \"\"\"Computes the min delays between all pairs of tasks in 'task_intervals'.\n",
    "\n",
    "    Args:\n",
    "      problem: The protobuf of the model.\n",
    "      task_intervals: The output of the AnalysePrecedenceGraph().\n",
    "\n",
    "    Returns:\n",
    "      a list of (task1, task2, min_delay_between_task1_and_task2)\n",
    "    \"\"\"\n",
    "    print(\"Computing the minimum delay between pairs of intervals\")\n",
    "    delays = {}\n",
    "    if (\n",
    "        problem.is_resource_investment\n",
    "        or problem.is_consumer_producer\n",
    "        or problem.is_rcpsp_max\n",
    "        or _DELAY_TIME_LIMIT.value <= 0.0\n",
    "    ):\n",
    "        return delays, None, False\n",
    "\n",
    "    complete_problem_assignment = None\n",
    "    num_optimal_delays = 0\n",
    "    num_delays_not_found = 0\n",
    "    optimal_found = True\n",
    "    for start_task, end_task, active_tasks in task_intervals:\n",
    "        min_delay, feasible_delay, assignment = SolveRcpsp(\n",
    "            problem,\n",
    "            \"\",\n",
    "            f\"num_search_workers:16,max_time_in_seconds:{_DELAY_TIME_LIMIT.value}\",\n",
    "            active_tasks,\n",
    "            start_task,\n",
    "            end_task,\n",
    "            [],\n",
    "            delays,\n",
    "        )\n",
    "        if min_delay != -1:\n",
    "            delays[(start_task, end_task)] = min_delay, feasible_delay\n",
    "            if start_task == 0 and end_task == len(problem.tasks) - 1:\n",
    "                complete_problem_assignment = assignment\n",
    "            if min_delay == feasible_delay:\n",
    "                num_optimal_delays += 1\n",
    "            else:\n",
    "                optimal_found = False\n",
    "        else:\n",
    "            num_delays_not_found += 1\n",
    "            optimal_found = False\n",
    "\n",
    "    print(f\"  - #optimal delays = {num_optimal_delays}\", flush=True)\n",
    "    if num_delays_not_found:\n",
    "        print(f\"  - #not computed delays = {num_delays_not_found}\", flush=True)\n",
    "\n",
    "    return delays, complete_problem_assignment, optimal_found\n",
    "\n",
    "\n",
    "def AcceptNewCandidate(problem, after, demand_map, current, candidate):\n",
    "    \"\"\"Check if candidate is compatible with the tasks in current.\"\"\"\n",
    "    for c in current:\n",
    "        if candidate in after[c] or c in after[candidate]:\n",
    "            return False\n",
    "\n",
    "    all_resources = range(len(problem.resources))\n",
    "    for res in all_resources:\n",
    "        resource = problem.resources[res]\n",
    "        if not resource.renewable:\n",
    "            continue\n",
    "        if (\n",
    "            sum(demand_map[(t, res)] for t in current) + demand_map[(candidate, res)]\n",
    "            > resource.max_capacity\n",
    "        ):\n",
    "            return False\n",
    "\n",
    "    return True\n",
    "\n",
    "\n",
    "def ComputePreemptiveLowerBound(problem, after, lower_bound):\n",
    "    \"\"\"Computes a preemtive lower bound for the makespan statically.\n",
    "\n",
    "    For this, it breaks all intervals into a set of intervals of size one.\n",
    "    Then it will try to assign all of them in a minimum number of configurations.\n",
    "    This is a standard complete set covering using column generation approach\n",
    "    where each column is a possible combination of itervals of size one.\n",
    "\n",
    "    Args:\n",
    "      problem: The probuf of the model.\n",
    "      after: a task to list of task dict that contains all tasks after a given\n",
    "        task.\n",
    "      lower_bound: A valid lower bound of the problem. It can be 0.\n",
    "\n",
    "    Returns:\n",
    "      a valid lower bound of the problem.\n",
    "    \"\"\"\n",
    "    # Check this is a single mode problem.\n",
    "    if (\n",
    "        problem.is_rcpsp_max\n",
    "        or problem.is_resource_investment\n",
    "        or problem.is_consumer_producer\n",
    "    ):\n",
    "        return lower_bound\n",
    "\n",
    "    demand_map = collections.defaultdict(int)\n",
    "    duration_map = {}\n",
    "    all_active_tasks = list(range(1, len(problem.tasks) - 1))\n",
    "    max_duration = 0\n",
    "    sum_of_demands = 0\n",
    "\n",
    "    for t in all_active_tasks:\n",
    "        task = problem.tasks[t]\n",
    "        if len(task.recipes) > 1:\n",
    "            return 0\n",
    "        recipe = task.recipes[0]\n",
    "        duration_map[t] = recipe.duration\n",
    "        for demand, resource in zip(recipe.demands, recipe.resources):\n",
    "            demand_map[(t, resource)] = demand\n",
    "            max_duration = max(max_duration, recipe.duration)\n",
    "            sum_of_demands += demand\n",
    "\n",
    "    print(\n",
    "        f\"Compute a bin-packing lower bound with {len(all_active_tasks)}\"\n",
    "        + \" active tasks\",\n",
    "        flush=True,\n",
    "    )\n",
    "    all_combinations = []\n",
    "\n",
    "    for t in all_active_tasks:\n",
    "        new_combinations = [[t]]\n",
    "\n",
    "        for c in all_combinations:\n",
    "            if AcceptNewCandidate(problem, after, demand_map, c, t):\n",
    "                new_combinations.append(c + [t])\n",
    "\n",
    "        all_combinations.extend(new_combinations)\n",
    "\n",
    "    print(f\"  - created {len(all_combinations)} combinations\")\n",
    "    if len(all_combinations) > 5000000:\n",
    "        return lower_bound  # Abort if too large.\n",
    "\n",
    "    # Solve the selection model.\n",
    "\n",
    "    # TODO(user): a few possible improvements:\n",
    "    # 1/  use \"dominating\" columns, i.e. if you can add a task to a column, then\n",
    "    #     do not use that column.\n",
    "    # 2/ Merge all task with exactly same demands into one.\n",
    "    model = cp_model.CpModel()\n",
    "    model.SetName(f\"lower_bound_{problem.name}\")\n",
    "\n",
    "    vars_per_task = collections.defaultdict(list)\n",
    "    all_vars = []\n",
    "    for c in all_combinations:\n",
    "        min_duration = max_duration\n",
    "        for t in c:\n",
    "            min_duration = min(min_duration, duration_map[t])\n",
    "        count = model.NewIntVar(0, min_duration, f\"count_{c}\")\n",
    "        all_vars.append(count)\n",
    "        for t in c:\n",
    "            vars_per_task[t].append(count)\n",
    "\n",
    "    # Each task must be performed.\n",
    "    for t in all_active_tasks:\n",
    "        model.Add(sum(vars_per_task[t]) >= duration_map[t])\n",
    "\n",
    "    # Objective\n",
    "    objective_var = model.NewIntVar(lower_bound, sum_of_demands, \"objective_var\")\n",
    "    model.Add(objective_var == sum(all_vars))\n",
    "\n",
    "    model.Minimize(objective_var)\n",
    "\n",
    "    # Solve model.\n",
    "    solver = cp_model.CpSolver()\n",
    "    solver.parameters.num_search_workers = 16\n",
    "    solver.parameters.max_time_in_seconds = _PREEMPTIVE_LB_TIME_LIMIT.value\n",
    "    status = solver.Solve(model)\n",
    "    if status == cp_model.OPTIMAL or status == cp_model.FEASIBLE:\n",
    "        status_str = \"optimal\" if status == cp_model.OPTIMAL else \"\"\n",
    "        lower_bound = max(lower_bound, int(solver.BestObjectiveBound()))\n",
    "        print(f\"  - {status_str} static lower bound = {lower_bound}\", flush=True)\n",
    "\n",
    "    return lower_bound\n",
    "\n",
    "\n",
    "def main(_):\n",
    "    rcpsp_parser = rcpsp.RcpspParser()\n",
    "    rcpsp_parser.parse_file(_INPUT.value)\n",
    "\n",
    "    problem = rcpsp_parser.problem()\n",
    "    PrintProblemStatistics(problem)\n",
    "\n",
    "    intervals_of_tasks, after = AnalyseDependencyGraph(problem)\n",
    "    delays, initial_solution, optimal_found = ComputeDelaysBetweenNodes(\n",
    "        problem, intervals_of_tasks\n",
    "    )\n",
    "\n",
    "    last_task = len(problem.tasks) - 1\n",
    "    key = (0, last_task)\n",
    "    lower_bound = delays[key][0] if key in delays else 0\n",
    "    if not optimal_found and _PREEMPTIVE_LB_TIME_LIMIT.value > 0.0:\n",
    "        lower_bound = ComputePreemptiveLowerBound(problem, after, lower_bound)\n",
    "\n",
    "    SolveRcpsp(\n",
    "        problem=problem,\n",
    "        proto_file=_OUTPUT_PROTO.value,\n",
    "        params=_PARAMS.value,\n",
    "        active_tasks=set(range(1, last_task)),\n",
    "        source=0,\n",
    "        sink=last_task,\n",
    "        intervals_of_tasks=intervals_of_tasks,\n",
    "        delays=delays,\n",
    "        in_main_solve=True,\n",
    "        initial_solution=initial_solution,\n",
    "        lower_bound=lower_bound,\n",
    "    )\n",
    "\n",
    "\n",
    "main()\n",
    "\n"
   ]
  }
 ],
 "metadata": {},
 "nbformat": 4,
 "nbformat_minor": 5
}
